!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!!
!! Simplified interface to PaStiX.
!!
!! \n
!!
!! This module provides a simplified interface to the <a
!! href="http://pastix.gforge.inria.fr/">PaStiX</a> solver, according
!! to the general layout given in \c mod_linsolver_base.
!!
!! The structure of this module is very similar to \c mod_mumpsintf.
!! In particular, each \c c_pastixpb is associated with one PaStiX
!! linear system; such linear system is built with a call to \c
!! factor, unless \c phase is present and equal to
!! <tt>"factorzation"</tt>, in which case it is assumed that the
!! system has already been created. The linear system is then
!! destroyed during a call to \c clean.
!!
!! The analysis phase corresponds to generating a new system and
!! setting the graph of this system. The factorization phase
!! corresponds to setting the matrix coefficients.
!!
!! \section parallel Parallel execution
!!
!! PaStiX has one centralized interface and two distributed
!! interfaces. All of them can be used for the parallel solution of a
!! linear system, the main difference being the input of the matrix:
!! centralized in one case, distributed in the other one. Here, we use
!! the <em>murge</em> distributed interface, since it is more general
!! and more similar to the one used for the other solvers. Its main
!! advantage is that it does not require that the matrix is
!! partitioned by columns among the processors <em>without
!! overlapping</em>, which would be a strong constraint.
!!
!! Relevant informations can be found <a
!! href="http://murge.gforge.inria.fr">here</a>, as well as in the
!! following files: <tt>example/src/fmurge.F90</tt>,
!! <tt>install/murge_int_double_real.in</tt>.
!!
!! \subsection repeated_entries Repeated entries
!!
!! The <em>murge</em> interface has two <em>modes</em>: \c
!! murge_assembly_fool and \c murge_assembly_respect. The second one
!! can be used when one is willing to respect the internal storage
!! scheme used by the solver, i.e., in the case of PaStiX, the column
!! partitioning without overlapping. This is of course the fastest
!! mode. The mode \c murge_assembly_fool gives more freedom, but
!! requires some communication overhead within the interface. Here, we
!! use this second mode. Summarizing, this amounts to:
!! <ul>
!!  <li> avoid building \c nodelist: this is essentially \c loc2glob
!!  which is used for the \c murge_assembly_respect mode; we don't use
!!  this here
!!  <li> always \c use murge_assembly_fool
!!  <li> use \c murge_assembly_add to add repeated contributions
!!  <li> use the following functions: murge_setrhs, murge_getsolution
!!  to deal with the right-hand-side locally without knowing the
!!  internal storage.
!! </ul>
!!
!! \section various Other comments
!!
!! The <em>murge</em> interface requires an initialization specifying
!! the "maximum number of solver instances that will be launched".
!! This is specified by a module parameter \c number_solver_instances.
module mod_pastixintf

!-----------------------------------------------------------------------

 use mod_utils, only: &
   t_realtime, my_second

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_sparse, only: &
   mod_sparse_initialized, &
   ! sparse types
   t_col, t_tri,&
   col2tri,     &
   transpose,   &
   clear

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_output_control, only: &
   mod_output_control_initialized, &
   elapsed_format, base_name

 use mod_linsolver_base, only: &
   mod_linsolver_base_initialized, &
   c_linpb

!-----------------------------------------------------------------------

 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_pastixintf_constructor, &
   mod_pastixintf_destructor,  &
   mod_pastixintf_initialized, &
   c_pastixpb

 private

!-----------------------------------------------------------------------

 include "murge.inc"

 ! Define some kind parameters for the murge interface
 integer(kind=murge_ints_kind) :: m_int
 real(kind=murge_coef_kind)    :: m_real
 integer, parameter :: mik = kind(m_int)
 integer, parameter :: mrk = kind(m_real)

 !> Linear PaStiX solver problem
 !!
 !! This type describes a linear system from the PaStiX viewpoint.
 type, extends(c_linpb), abstract :: c_pastixpb
  !> solve \f$A^Tx=b\f$
  logical :: transposed_mat = .false.
  !> matrix size
  integer :: gn
  !> local matrix, with zero based indexing (as always with \c t_col)
  type(t_col), pointer :: m
  !> local to global map, 0-based (as in \c mod_sparse)
  !!
  !! \note The lower bound must be zero, and the entries of the array
  !! must be zero-based indexes in the global matrix.
  integer, pointer :: gij(:)
  !> local right-hand side
  real(wp), pointer :: rhs(:)
  !> MPI communicator
  integer :: mpi_comm
  !> size of \c gij
  integer(mik), private :: nloc
  !> copy of \c gij with the proper kind type parameter
  integer(mik), private, allocatable :: c_gij(:)
  !> solver instance
  integer(mik), private :: sol_ins
  logical, private :: sys_set   = .false. !< internal consistency check
  logical, private :: coeff_set = .false. !< used in \c pastix_factor
  logical, private :: rhs_set   = .false. !< used in \c pastix_solve
 contains
  procedure, pass(s) :: factor => pastix_factor
  procedure, pass(s) :: solve  => pastix_solve
  procedure, pass(s) :: clean  => pastix_clean
  procedure, nopass :: working_implementation => pastix_wi
  procedure(i_xassign), deferred, pass(s) :: xassign
 end type c_pastixpb

 !> Convert the PaStiX solution into a \c c_stv object
 abstract interface
  subroutine i_xassign(x,s,pastix_x)
   import :: wp, c_stv, c_pastixpb
   implicit none
   real(wp),         intent(in) :: pastix_x(:)
   class(c_pastixpb), intent(inout) :: s
   class(c_stv),     intent(inout) :: x
  end subroutine i_xassign
 end interface

 !> Maximum number of solver instances
 integer(mik), parameter :: number_solver_instances = 5
 !> Counter of the active solver instances
 integer(mik) :: active_sol_inst

 logical, protected ::               &
   mod_pastixintf_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_pastixintf'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_pastixintf_constructor()

  integer(mik) :: ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
          (mod_kinds_initialized.eqv..false.) .or. &
         (mod_sparse_initialized.eqv..false.) .or. &
      (mod_mpi_utils_initialized.eqv..false.) .or. &
 (mod_output_control_initialized.eqv..false.) .or. &
 (mod_linsolver_base_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_pastixintf_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Starting MURGE
   call murge_initialize( number_solver_instances , ierr )
   if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
     'Problems initializing the murge interface')

   ! initialize the solver instance counter (numbered from 0)
   active_sol_inst = -1

   mod_pastixintf_initialized = .true.
 end subroutine mod_pastixintf_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_pastixintf_destructor()

  integer(mik) :: ierr
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_pastixintf_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   call murge_finalize(ierr)
   if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
     'Problems finalizing the murge interface')

   mod_pastixintf_initialized = .false.
 end subroutine mod_pastixintf_destructor

!-----------------------------------------------------------------------

  pure function pastix_wi()
   logical :: pastix_wi
    pastix_wi = .true.
  end function pastix_wi

!-----------------------------------------------------------------------

 !> Factor matrix \c m
 !!
 !! The easiest way to insert the matrix elements with the
 !! <em>murge</em> interface is converting the matrix to \c t_tri and
 !! inserting the elements one-by-one.
 !!
 !! \note In the numerical factorization step, during the insertion of
 !! the matrix coefficients, we use \c murge_assembly_add to add
 !! repeated entries. This is fine if the matrix is being created from
 !! scratch, but is the matrix is being redefined one must first set
 !! the pre-existing values to zero.
 subroutine pastix_factor(s,phase)
  class(c_pastixpb), intent(inout) :: s
  character(len=*), intent(in), optional :: phase
 
  logical :: do_analysis, do_factorization
  integer :: i
  integer(mik) :: ierr
  type(t_tri) :: m_tri
  character(len=*), parameter :: &
    this_sub_name = 'pastix_factor'

   !--------------------------------------------------------------------
   ! 0) Preliminaries

   ! define the required action
   do_analysis = .true.; do_factorization = .true. ! default: both
   if(present(phase)) then
     select case(trim(phase))
      case('analysis')
       do_factorization = .false.
      case('factorization')
       do_analysis      = .false.
      case default
       call error(this_sub_name,this_mod_name, &
           'Unknown phase "'//trim(phase)//'"' )
     end select
   endif

   ! convert the matrix to tri format
   if(s%transposed_mat) then
     m_tri = transpose(col2tri(s%m))
   else
     m_tri =           col2tri(s%m)
   endif

   !--------------------------------------------------------------------
   ! 1) Analysis phase

   analysis_if: if(do_analysis) then

     if(s%sys_set) call error(this_sub_name,this_mod_name,           &
       'An analysis exists already for this system: call "clean" '// &
       'before initializing a new linear problem.')

     !------------------------------------------------------------------
     ! 1.1) Solver setup

     ! create a new solver instance
     active_sol_inst = active_sol_inst + 1
     s%sol_ins = active_sol_inst

     ! set the default options
     call murge_setdefaultoptions(s%sol_ins , 0_mik , ierr)

     ! generic, unsymmetric matrix
     call murge_setoptionint(s%sol_ins, murge_iparam_sym, &
                             murge_boolean_false, ierr)

     ! use zero-based indexes
     call murge_setoptionint(s%sol_ins,murge_iparam_baseval,0_mik,ierr)

     !------------------------------------------------------------------
     ! 1.2) Graph setup

     s%nloc = size(s%gij)
     allocate(s%c_gij(0:s%nloc-1))
     s%c_gij = int( s%gij , mik ) ! this work even for lbound(gij).ne.0
   
     call murge_graphbegin(s%sol_ins, s%gn, m_tri%nz, ierr)
     if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
       'Problems calling murge_graphbegin')

     do i=1,m_tri%nz
       call murge_graphedge(s%sol_ins , &
              row = s%c_gij( m_tri%ti(i) ) , &
              col = s%c_gij( m_tri%tj(i) ) , &
           ierror = ierr)
       if(ierr.ne.murge_success)call error(this_sub_name,this_mod_name,&
         'Problems inserting a graph-edge with murge_graphedge')
     enddo

     ! close the graph
     call murge_graphend(s%sol_ins , ierr)
     if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
       'Problems calling murge_graphend')

     s%sys_set = .true.
   endif analysis_if

   !--------------------------------------------------------------------
   ! 2) Factorization phase

   factorization_if: if(do_factorization) then

     if(.not.s%sys_set) call error(this_sub_name,this_mod_name,      &
       'An analysis must be performed before doing the numerical '// &
       'factorizazion.')

     !------------------------------------------------------------------
     ! 2.0) Reset a pre-existing matrix if necessary

     if(s%coeff_set) then
       call murge_assemblybegin(s%sol_ins, s%gn, m_tri%nz,             &
           murge_assembly_ovw, murge_assembly_ovw, murge_assembly_fool,&
           murge_boolean_false, ierr)

       do i=1,m_tri%nz
         call murge_assemblysetvalue(s%sol_ins , &
                row = s%c_gij( m_tri%ti(i) ) ,   &
                col = s%c_gij( m_tri%tj(i) ) ,   &
              value = 0.0_mrk,                   &
             ierror = ierr)
         if(ierr.ne.murge_success) call error(this_sub_name, &
            this_mod_name, 'Problems setting a value to zero')
       enddo
     endif

     !------------------------------------------------------------------
     ! 2.1) Matrix assembly
  
     call murge_assemblybegin(s%sol_ins, s%gn, m_tri%nz,              &
         murge_assembly_add, murge_assembly_add, murge_assembly_fool, &
         murge_boolean_false, ierr)
  
     do i=1,m_tri%nz
       call murge_assemblysetvalue(s%sol_ins , &
              row = s%c_gij( m_tri%ti(i) )   , &
              col = s%c_gij( m_tri%tj(i) )   , &
            value = real(  m_tri%tx(i) ,mrk) , &
           ierror = ierr)
       if(ierr.ne.murge_success)call error(this_sub_name,this_mod_name,&
         'Problems inserting a value')
     enddo
  
     call murge_assemblyend(s%sol_ins, ierr)
     if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
       'Problems calling murge_assemblyend')

     s%coeff_set = .true.

   endif factorization_if
  
   !--------------------------------------------------------------------
   ! 3) Finalization
   call clear(m_tri)

 end subroutine pastix_factor

!-----------------------------------------------------------------------

 subroutine pastix_solve(x,s)
  class(c_pastixpb), intent(inout) :: s
  class(c_stv),     intent(inout) :: x

  integer(mik) :: ierr
  real(mrk), allocatable :: x_array(:)
  character(len=*), parameter :: &
    this_sub_name = 'pastix_solve'
 
   if(.not.s%sys_set) call error(this_sub_name,this_mod_name,   &
     'Trying to solve a system while the matrix is not defined.'   )
   if(.not.s%coeff_set) call error(this_sub_name,this_mod_name, &
     'Trying to solve a system while the matrix is not factorized.')

   allocate(x_array(s%nloc))

   !--------------------------------------------------------------------
   ! 1) Set the RHS

   ! Zeroth previous values is required
   if(s%rhs_set) then
     x_array = 0.0_mrk
     call murge_setrhs(s%sol_ins, &
            n = s%nloc , & ! number of coeffs to set
     coefsidx = s%c_gij, & ! coeff list
            b = x_array, &
           op = murge_assembly_ovw, &
          op2 = murge_assembly_ovw, &
         mode = murge_assembly_fool, ierror = ierr)
     if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
       'Problems calling murge_setrhs to zeroth previous values.')
   endif

   call murge_setrhs(s%sol_ins, &
            n =       s%nloc ,      & ! number of coeffs to set
     coefsidx =       s%c_gij,      & ! coeff list
            b = real( s%rhs  ,mrk), & ! rhs entries
           op = murge_assembly_add, &
          op2 = murge_assembly_add, &
         mode = murge_assembly_fool, ierror = ierr)
   if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
     'Problems calling murge_setrhs to set the rhs.')

   s%rhs_set = .true.

   !--------------------------------------------------------------------
   ! 2) Solve the system

   call murge_getsolution(s%sol_ins, &
            n = s%nloc  , & ! number of coeffs to set
     coefsidx = s%c_gij , & ! coeff list
            x = x_array , & ! solution
         mode = murge_assembly_fool, ierror = ierr)
   if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
     'Problems calling murge_getsolution.')

   !--------------------------------------------------------------------
   ! 3) Copy back the solution
   call s%xassign(x,real(x_array,wp))

   deallocate(x_array)
 end subroutine pastix_solve

!-----------------------------------------------------------------------
 
 subroutine pastix_clean(s)
  class(c_pastixpb), intent(inout) :: s

  integer(mik) :: ierr
  character(len=*), parameter :: &
    this_sub_name = 'pastix_clean'

   if(s%sys_set) then

     call murge_clean(s%sol_ins, ierr)
     if(ierr.ne.murge_success) call error(this_sub_name,this_mod_name, &
       'Problems calling murge_clean.')

     deallocate(s%c_gij)
 
     ! decrement the instance counter
     active_sol_inst = active_sol_inst - 1

   endif

   s%sys_set = .false.

 end subroutine pastix_clean

!-----------------------------------------------------------------------
 
end module mod_pastixintf

